library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.3.2
## Warning: package 'tidyr' was built under R version 4.3.2
library(here)
library(phyloseq)
library(vegan)
## Warning: package 'lattice' was built under R version 4.3.2
theme_set(theme_bw())
filt_fun <- function (x, min_reads = 100, min_samples = 5) {
(sum(x) > min_reads) & (sum(x > 0) > min_samples)
}
ps <- readRDS(here('data','following_study','ps_rarefied.rds')) %>%
filter_taxa(filt_fun, prune = TRUE)
sample_data(ps)[is.na(sample_data(ps))] <- 'NA'
sam <- data.frame(sample_data(ps))
max.core <- parallel::detectCores()
alpha.diversity <- data.frame(sample_data(ps),estimate_richness(ps),read_depth = rowSums(otu_table(ps)))
ggplot(alpha.diversity,aes(x = Household, y = Shannon, group = Household)) +
geom_point() + geom_line()
anova(lm(Shannon~Household, data = alpha.diversity))
## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Household 48 8.4359 0.17575 1.566 0.06064 .
## Residuals 49 5.4992 0.11223
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(alpha.diversity,aes(x = Epileptic.or.Control, y = Shannon)) +
geom_boxplot() + geom_jitter(height = 0, width = 0.25)
Paired t test
epileptic <- alpha.diversity %>% filter(Epileptic.or.Control == 'Epileptic')
control <- alpha.diversity %>% filter(Epileptic.or.Control == 'Control')
t.test(epileptic$Shannon, control$Shannon, paired = TRUE)
##
## Paired t-test
##
## data: epileptic$Shannon and control$Shannon
## t = -0.22662, df = 48, p-value = 0.8217
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
## -0.1529073 0.1219302
## sample estimates:
## mean difference
## -0.01548853
ANOVA with Household effect added
anova(lm(Shannon~Household + Epileptic.or.Control, data = alpha.diversity))
## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Household 48 8.4359 0.175747 1.5357 0.07042 .
## Epileptic.or.Control 1 0.0059 0.005877 0.0514 0.82168
## Residuals 48 5.4933 0.114444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(alpha.diversity) +
geom_point(aes(x = Breed.Group..1., y = Shannon, colour = Breed.Group..1.)) +
facet_wrap(~Epileptic.or.Control) +
theme(axis.text.x = element_blank(), axis.ticks.x.bottom = element_blank())
# here breed group with na is removed
anova(lm(Shannon~Household + Breed.Group..1., data = alpha.diversity %>% filter(Breed.Group..1. != 'NA')))
## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Household 45 7.2456 0.16101 1.6626 0.05393 .
## Breed.Group..1. 4 1.3804 0.34511 3.5637 0.01431 *
## Residuals 39 3.7768 0.09684
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
alpha.diversity.epi <- alpha.diversity %>% filter(Epileptic.or.Control == 'Epileptic')
ggplot(alpha.diversity.epi, aes(x = Pheno.Y.N, y = Shannon)) +
geom_boxplot() + geom_jitter(height = 0, width = 0.25)
t.test(Shannon~Pheno.Y.N, data = alpha.diversity.epi)
##
## Welch Two Sample t-test
##
## data: Shannon by Pheno.Y.N
## t = 0.84674, df = 33.509, p-value = 0.4031
## alternative hypothesis: true difference in means between group No and group Yes is not equal to 0
## 95 percent confidence interval:
## -0.1450834 0.3521435
## sample estimates:
## mean in group No mean in group Yes
## 2.730402 2.626872
ggplot(alpha.diversity, aes(x = Sex, y = Shannon)) +
geom_boxplot() + geom_jitter(height = 0, width = 0.25)
anova(lm(Shannon~Household + Sex, data = alpha.diversity))
## Analysis of Variance Table
##
## Response: Shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Household 48 8.4359 0.17575 1.6253 0.04789 *
## Sex 1 0.3088 0.30882 2.8560 0.09752 .
## Residuals 48 5.1903 0.10813
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
t.test(Shannon~Sex, data = alpha.diversity)
##
## Welch Two Sample t-test
##
## data: Shannon by Sex
## t = -1.315, df = 87.559, p-value = 0.1919
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
## -0.25315537 0.05154705
## sample estimates:
## mean in group F mean in group M
## 2.658970 2.759775
ggplot(alpha.diversity, aes(x = as.numeric(Age..months.), y = Shannon)) +
geom_point()
## Warning in FUN(X[[i]], ...): NAs introduced by coercion
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
lm(Shannon~as.numeric(Age..months.), data = alpha.diversity) %>% summary()
## Warning in eval(predvars, data, env): NAs introduced by coercion
##
## Call:
## lm(formula = Shannon ~ as.numeric(Age..months.), data = alpha.diversity)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.92260 -0.28131 -0.01761 0.27021 0.78472
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.6317958 0.0896466 29.357 <2e-16 ***
## as.numeric(Age..months.) 0.0007684 0.0010859 0.708 0.481
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3742 on 94 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.005298, Adjusted R-squared: -0.005284
## F-statistic: 0.5007 on 1 and 94 DF, p-value: 0.4809
ord <- ps %>% ordinate(method = "MDS", distance = 'unifrac')
plot_ordination(ps, ord)
The PCoA plot shows a very separate pattern of four clusters. Let’s take a look at the phlogenetic tree of our data.
plot_tree(ps, "treeonly")
Some of the taxa have extremely long branch. Are they influential?
# Warning: Not all nodes are descendants in this tree
if (require(adephylo)) {dist.root <- adephylo::distRoot(phy_tree(ps))}
## Loading required package: adephylo
## Loading required package: ade4
## Found more than one class "phylo" in cache; using the first, from namespace 'phyloseq'
## Also defined by 'RNeXML'
hist(dist.root); dist.root %>% sort(decreasing = TRUE) %>% head(20)
## ASV601 ASV438 ASV563 ASV287 ASV484 ASV133 ASV562 ASV437
## 20.819315 5.082379 4.338564 4.267155 4.245708 4.159508 4.154919 4.121791
## ASV450 ASV286 ASV478 ASV357 ASV446 ASV454 ASV423 ASV380
## 4.088335 4.077083 4.074518 4.060563 4.045067 4.041199 4.026320 1.631610
## ASV213 ASV644 ASV268 ASV41
## 1.617836 1.569113 1.417223 1.411706
detach('package:adephylo', character.only = TRUE, unload = TRUE, force = TRUE)
Let’s see which ASV has distance from root longer than 10.
long.branch <- dist.root[dist.root > 5]
long.branch
## ASV601 ASV438
## 20.819315 5.082379
Are these ASVs special?
dada2:::pfasta(refseq(ps)[names(long.branch)], ids = names(long.branch))
## >ASV601
## GGAGAGCCATGTCTTATAAGTATACCATACTCTTTTAACTTCTTAATGTTCTTTTCTATGCTTCGTATTCGAAAGAAATAAATTCACATTTTTTTCTGTGTATCAGTTAATTCATAATCCTGAAATCTCGGCTTCGATAGTCCATACTCTCTCTTATCCAGAAACACTGCACGATCTGTTCCCTTGAATACCGCACATTGTGTTTTAGCGAATGAAAAATAGTCCGATGTCAGCAACGCATAGGCATTCGTTGCCAGAAGCTGTCCCTCACTCTGCTTTAAAATCTTCCACCAACATC
## >ASV438
## CAAATCCGAACTGCAAAATACTCATGGAAAGAATCGCTCCTGTTTTTATTCCTTCTCCTAACATGGTAAGAACAGAAACTGGCACTCCAAATACAGAAAAACCACAAAGCATACAAATAAAAAGCCATCCGCCTCTCTGTGCCAGAACCTGGAAAAAATATTCTTTTCCAGAAACATCTTCCCCCGCAAAAGCACCCAAAAGGTAAAATACATGGACTGTTTCCTGTTTCCACTGCATCTTCCACAAAAGATTTGGAAGAACCGTTCCC
Result returned from BLAST (web interface):
| ASV | Taxonomy | % identity |
|---|---|---|
| ASV601 | Anaerobutyricum hallii | 95.1 |
| ASV438 | uncultured Blautia sp. | 97.3 |
Now we remove this outlier ASVs
all.taxa <- taxa_names(ps)
keep.taxa <- all.taxa[!(all.taxa %in% names(long.branch))]
ps.sub <- ps %>% prune_taxa(keep.taxa, .)
Convert ps.sub to proportion
ps.sub <- ps.sub %>% transform_sample_counts(function(x) x / sum(x))
sample_sums(ps.sub) %>% summary()
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 1 1 1 1 1
plot_ord <- function(data, factor, method, distance){
data.ord <- ordinate(data, method = method, distance = distance)
p <- plot_ordination(data, data.ord, color = factor)
p <- p + stat_ellipse(type = "t",geom = "polygon",alpha = 0)
p <- p + ggtitle(str_c(factor,method,distance, sep = ' - '))
print(p)
}
permanova <- function(data, formula, method, permutations=1e4, strata = NULL, core = max.core){
dist.matrix <- phyloseq::distance(data, method=method)
df <- data.frame(sample_data(data))
model <- as.formula(paste0('dist.matrix~', formula))
if (!is_null(strata)) {strata <- df[,strata]}
result <- adonis2(model,
data = df,
permutations=permutations,
strata = strata,
parallel = core)
return(result)
}
permdisp <- function(data, group, method, permutations=1e4, pairwise = FALSE, core = max.core){
dist.matrix <- phyloseq::distance(data, method=method)
df <- data.frame(sample_data(data))
beta.disp <- betadisper(dist.matrix, group = df[,group])
result <- permutest(beta.disp, permutations = permutations, pairwise = pairwise, type = 'centroid')
return(result)
}
permanova(ps.sub, 'Household', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 11.7188 0.68925 2.2643 9.999e-05 ***
## Residual 49 5.2834 0.31075
## Total 97 17.0021 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot_ord(ps.sub, 'Epileptic.or.Control','MDS','bray')
plot_ord(ps.sub, 'Epileptic.or.Control','MDS','wunifrac')
plot_ord(ps.sub, 'Epileptic.or.Control','NMDS','bray')
## Run 0 stress 0.2084864
## Run 1 stress 0.2041351
## ... New best solution
## ... Procrustes: rmse 0.05361631 max resid 0.4049279
## Run 2 stress 0.2035787
## ... New best solution
## ... Procrustes: rmse 0.03136479 max resid 0.2746453
## Run 3 stress 0.2029836
## ... New best solution
## ... Procrustes: rmse 0.01769698 max resid 0.07861787
## Run 4 stress 0.208685
## Run 5 stress 0.2030752
## ... Procrustes: rmse 0.004255481 max resid 0.02643523
## Run 6 stress 0.2088515
## Run 7 stress 0.2064062
## Run 8 stress 0.2162702
## Run 9 stress 0.2029817
## ... New best solution
## ... Procrustes: rmse 0.0006812619 max resid 0.004905256
## ... Similar to previous best
## Run 10 stress 0.2040754
## Run 11 stress 0.2100049
## Run 12 stress 0.2053367
## Run 13 stress 0.2113234
## Run 14 stress 0.2106682
## Run 15 stress 0.2312135
## Run 16 stress 0.2039773
## Run 17 stress 0.2041273
## Run 18 stress 0.2122018
## Run 19 stress 0.209445
## Run 20 stress 0.2047825
## *** Best solution repeated 1 times
plot_ord(ps.sub, 'Epileptic.or.Control','NMDS','wunifrac')
## Run 0 stress 0.1726183
## Run 1 stress 0.1740023
## Run 2 stress 0.1762893
## Run 3 stress 0.1718306
## ... New best solution
## ... Procrustes: rmse 0.0235026 max resid 0.173255
## Run 4 stress 0.1872898
## Run 5 stress 0.1706236
## ... New best solution
## ... Procrustes: rmse 0.01010811 max resid 0.09705579
## Run 6 stress 0.1718317
## Run 7 stress 0.1861904
## Run 8 stress 0.1899837
## Run 9 stress 0.1706237
## ... Procrustes: rmse 5.413324e-05 max resid 0.0003045123
## ... Similar to previous best
## Run 10 stress 0.1726201
## Run 11 stress 0.1696998
## ... New best solution
## ... Procrustes: rmse 0.02436172 max resid 0.1715287
## Run 12 stress 0.1947747
## Run 13 stress 0.1822626
## Run 14 stress 0.1876531
## Run 15 stress 0.1802945
## Run 16 stress 0.1864107
## Run 17 stress 0.1865221
## Run 18 stress 0.1856708
## Run 19 stress 0.1891249
## Run 20 stress 0.185231
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 1: no. of iterations >= maxit
## 19: stress ratio > sratmax
permanova(ps.sub, 'Epileptic.or.Control', 'bray', strata = 'Household')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Epileptic.or.Control 1 0.1499 0.00882 0.854 0.111
## Residual 96 16.8522 0.99118
## Total 97 17.0021 1.00000
permanova(ps.sub, 'Epileptic.or.Control', 'wunifrac', strata = 'Household')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Blocks: strata
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Epileptic.or.Control 1 0.1149 0.01615 1.5756 0.0189 *
## Residual 96 6.9986 0.98385
## Total 97 7.1135 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permdisp(ps.sub, 'Epileptic.or.Control', 'bray')
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.00251 0.0025142 0.1605 10000 0.6924
## Residuals 96 1.50379 0.0156645
permdisp(ps.sub, 'Epileptic.or.Control', 'wunifrac')
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.02183 0.0218264 2.2935 10000 0.13
## Residuals 96 0.91358 0.0095165
a table of # of breed. frequency table
(data.frame(sample_data(ps.sub))$Breed.Group..1.) %>% table()
## .
## Herder NA Pointer Spaniel Retriever Scent Hound
## 20 9 26 26 5
## Sight Hound Sled Dog Terrier
## 3 6 3
plot_ord(ps.sub, 'Breed.Group..1.','MDS','bray')
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
plot_ord(ps.sub, 'Breed.Group..1.','MDS','wunifrac')
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
plot_ord(ps.sub, 'Breed.Group..1.','NMDS','bray')
## Run 0 stress 0.2084864
## Run 1 stress 0.2040108
## ... New best solution
## ... Procrustes: rmse 0.05685435 max resid 0.4005662
## Run 2 stress 0.2089403
## Run 3 stress 0.2093198
## Run 4 stress 0.20545
## Run 5 stress 0.2173889
## Run 6 stress 0.219423
## Run 7 stress 0.203331
## ... New best solution
## ... Procrustes: rmse 0.03779231 max resid 0.3175873
## Run 8 stress 0.2047923
## Run 9 stress 0.2042308
## Run 10 stress 0.2047717
## Run 11 stress 0.2049309
## Run 12 stress 0.2037882
## ... Procrustes: rmse 0.03544949 max resid 0.3059876
## Run 13 stress 0.2047891
## Run 14 stress 0.2049428
## Run 15 stress 0.2041995
## Run 16 stress 0.2030192
## ... New best solution
## ... Procrustes: rmse 0.01241197 max resid 0.0799837
## Run 17 stress 0.2041269
## Run 18 stress 0.2029852
## ... New best solution
## ... Procrustes: rmse 0.005615208 max resid 0.03165945
## Run 19 stress 0.211464
## Run 20 stress 0.2040137
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 5: no. of iterations >= maxit
## 15: stress ratio > sratmax
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
plot_ord(ps.sub, 'Breed.Group..1.','NMDS','wunifrac')
## Run 0 stress 0.1726183
## Run 1 stress 0.184764
## Run 2 stress 0.1697108
## ... New best solution
## ... Procrustes: rmse 0.01847867 max resid 0.1443311
## Run 3 stress 0.1747189
## Run 4 stress 0.1724669
## Run 5 stress 0.1726183
## Run 6 stress 0.1724713
## Run 7 stress 0.1872728
## Run 8 stress 0.1706105
## Run 9 stress 0.1921029
## Run 10 stress 0.1791179
## Run 11 stress 0.1820495
## Run 12 stress 0.1853757
## Run 13 stress 0.1800669
## Run 14 stress 0.1921885
## Run 15 stress 0.1849801
## Run 16 stress 0.1773687
## Run 17 stress 0.1884582
## Run 18 stress 0.1819956
## Run 19 stress 0.1792786
## Run 20 stress 0.1912135
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 19: stress ratio > sratmax
## 1: scale factor of the gradient < sfgrmin
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
# same breed group v.s. diff breed group in household
# check if data are paired within household
if (!identical(sample_data(ps)$Household[seq(1,98,2)],
sample_data(ps)$Household[seq(1,98,2)+1])) stop('data is not paired')
same.breed.group <- table(sample_data(ps)$Household, sample_data(ps)$Breed.Group..1.) %>%
apply(1, function(x) any(x == 2)) # if dog's grouping are same with in each house
same.breed.group;sum(same.breed.group)
## 1 10 11 12 13 14 15 16 17 18 19 2 20
## FALSE TRUE TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE TRUE TRUE
## 21 22 23 24 25 26 27 28 29 3 30 31 32
## TRUE FALSE TRUE FALSE TRUE TRUE FALSE TRUE FALSE TRUE TRUE TRUE TRUE
## 33 34 35 36 38 39 4 40 41 42 43 44 45
## TRUE TRUE TRUE TRUE FALSE FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## 47 48 49 5 50 51 6 7 8 9
## TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE FALSE
## [1] 39
Here we see 39 out of 49 household have dogs in same breed group
permanova(ps.sub, 'Household + Breed.Group..1.', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 11.7188 0.68925 2.3509 9.999e-05 ***
## Breed.Group..1. 6 0.8179 0.04811 1.3127 0.06899 .
## Residual 43 4.4655 0.26264
## Total 97 17.0021 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.sub, 'Household + Breed.Group..1.', 'wunifrac')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 4.8365 0.67991 2.2564 9.999e-05 ***
## Breed.Group..1. 6 0.3568 0.05015 1.3316 0.1275
## Residual 43 1.9202 0.26993
## Total 97 7.1135 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Here we
ps.control <- ps.sub %>% subset_samples(Epileptic.or.Control == 'Control')
ps.epileptic <- ps.sub %>% subset_samples(Epileptic.or.Control == 'Epileptic')
perm.breed.bray.ctl <- permanova(ps.control, 'Breed.Group..1.', 'bray')
perm.breed.wunif.ctl <- permanova(ps.control, 'Breed.Group..1.', 'wunifrac')
perm.breed.bray.ctl;perm.breed.wunif.ctl
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Breed.Group..1. 7 1.6962 0.19479 1.4169 0.0319 *
## Residual 41 7.0114 0.80521
## Total 48 8.7076 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Breed.Group..1. 7 0.8423 0.21377 1.5925 0.032 *
## Residual 41 3.0977 0.78623
## Total 48 3.9400 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
perm.breed.bray.epi <- permanova(ps.epileptic, 'Breed.Group..1.', 'bray')
perm.breed.wunif.epi <- permanova(ps.epileptic, 'Breed.Group..1.', 'wunifrac')
perm.breed.bray.epi;perm.breed.wunif.epi
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Breed.Group..1. 7 1.2777 0.15688 1.0898 0.2615
## Residual 41 6.8669 0.84312
## Total 48 8.1446 1.00000
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Breed.Group..1. 7 0.55214 0.18052 1.2902 0.1023
## Residual 41 2.50647 0.81948
## Total 48 3.05861 1.00000
\[X^2=-2 \sum_{i=1}^k \ln \left(p_i\right)\]
bray.p <- c(perm.breed.bray.ctl$`Pr(>F)`[1], perm.breed.bray.epi$`Pr(>F)`[1])
# Calculate the combined p-value with fisher's method
pchisq( -2 * sum(log(bray.p)), df = 4, lower.tail = FALSE)
## [1] 0.04826188
wunif.p <- c(perm.breed.wunif.ctl$`Pr(>F)`[1], perm.breed.wunif.epi$`Pr(>F)`[1])
# Calculate the combined p-value with fisher's method
pchisq( -2 * sum(log(wunif.p)), df = 4, lower.tail = FALSE)
## [1] 0.02200095
With using the weighted-unifrac distance, there’s at
least one of the null hypothesis is rejected.
Here we are only focusing on epileptic dogs.
plot_ord(ps.epileptic, 'Pheno.Y.N','MDS','bray')
plot_ord(ps.epileptic, 'Pheno.Y.N','MDS','wunifrac')
plot_ord(ps.epileptic, 'Pheno.Y.N','NMDS','bray')
## Run 0 stress 0.2030888
## Run 1 stress 0.2049861
## Run 2 stress 0.2059649
## Run 3 stress 0.2101223
## Run 4 stress 0.2165297
## Run 5 stress 0.1990618
## ... New best solution
## ... Procrustes: rmse 0.0459343 max resid 0.1625879
## Run 6 stress 0.1999106
## Run 7 stress 0.2132579
## Run 8 stress 0.2146445
## Run 9 stress 0.2087338
## Run 10 stress 0.2049879
## Run 11 stress 0.2020433
## Run 12 stress 0.2239229
## Run 13 stress 0.2080408
## Run 14 stress 0.2194199
## Run 15 stress 0.2106283
## Run 16 stress 0.2055968
## Run 17 stress 0.2177038
## Run 18 stress 0.2281718
## Run 19 stress 0.206055
## Run 20 stress 0.2180574
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
plot_ord(ps.epileptic, 'Pheno.Y.N','NMDS','wunifrac')
## Run 0 stress 0.1857
## Run 1 stress 0.1942337
## Run 2 stress 0.1854942
## ... New best solution
## ... Procrustes: rmse 0.07777616 max resid 0.237673
## Run 3 stress 0.1873453
## Run 4 stress 0.1874356
## Run 5 stress 0.19423
## Run 6 stress 0.1906979
## Run 7 stress 0.1927119
## Run 8 stress 0.1814438
## ... New best solution
## ... Procrustes: rmse 0.0508053 max resid 0.2199367
## Run 9 stress 0.1831429
## Run 10 stress 0.1918678
## Run 11 stress 0.1932858
## Run 12 stress 0.4028392
## Run 13 stress 0.193139
## Run 14 stress 0.1884053
## Run 15 stress 0.1857624
## Run 16 stress 0.1943164
## Run 17 stress 0.1937669
## Run 18 stress 0.2090482
## Run 19 stress 0.1830395
## Run 20 stress 0.1894807
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
permanova(ps.epileptic, 'Pheno.Y.N','bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Pheno.Y.N 1 0.2090 0.02566 1.2379 0.1963
## Residual 47 7.9356 0.97434
## Total 48 8.1446 1.00000
permanova(ps.epileptic, 'Pheno.Y.N','wunifrac')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Pheno.Y.N 1 0.06486 0.02121 1.0182 0.4158
## Residual 47 2.99375 0.97879
## Total 48 3.05861 1.00000
permdisp(ps.epileptic, 'Pheno.Y.N','bray')
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.01951 0.019507 1.437 10000 0.2338
## Residuals 47 0.63801 0.013575
permdisp(ps.epileptic, 'Pheno.Y.N','wunifrac')
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.03162 0.0316221 4.2019 10000 0.0467 *
## Residuals 47 0.35370 0.0075256
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# within epileptic dogs
table(sample_data(ps)$Seizure.Freedom..Y.N., sample_data(ps)$Epileptic.or.Control)
##
## Control Epileptic
## NA 49 3
## No 0 32
## Yes 0 14
table(sample_data(ps)$Seizure.Control..Satisfactory.Unsatisfactory., sample_data(ps)$Epileptic.or.Control)
##
## Control Epileptic
## NA 49 3
## S 0 33
## US 0 13
Does health condition related to dog’s sex? Here we conduct a Chi-squared test.
# chisq test on epi male dog v.s. female.
# not within health condition
tb <- table(sample_data(ps)$Sex, sample_data(ps)$Epileptic.or.Control)
tb
##
## Control Epileptic
## F 33 25
## M 16 24
chisq.test(tb)
##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: tb
## X-squared = 2.0698, df = 1, p-value = 0.1502
plot_ord(ps, 'Sex','MDS','bray')
plot_ord(ps, 'Sex','MDS','wunifrac')
plot_ord(ps, 'Sex','NMDS','bray')
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.2261393
## Run 1 stress 0.2450121
## Run 2 stress 0.2605166
## Run 3 stress 0.2507737
## Run 4 stress 0.2445078
## Run 5 stress 0.2364009
## Run 6 stress 0.2282095
## Run 7 stress 0.229298
## Run 8 stress 0.2365472
## Run 9 stress 0.2371477
## Run 10 stress 0.2342887
## Run 11 stress 0.2296419
## Run 12 stress 0.2547756
## Run 13 stress 0.2597614
## Run 14 stress 0.2361598
## Run 15 stress 0.2432431
## Run 16 stress 0.2491273
## Run 17 stress 0.235258
## Run 18 stress 0.2485772
## Run 19 stress 0.2446727
## Run 20 stress 0.2526808
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress ratio > sratmax
plot_ord(ps, 'Sex','NMDS','wunifrac')
## Run 0 stress 0.1724074
## Run 1 stress 0.1810284
## Run 2 stress 0.1705573
## ... New best solution
## ... Procrustes: rmse 0.0224693 max resid 0.1611144
## Run 3 stress 0.1808744
## Run 4 stress 0.1971364
## Run 5 stress 0.189983
## Run 6 stress 0.1817739
## Run 7 stress 0.1836256
## Run 8 stress 0.1919375
## Run 9 stress 0.1854364
## Run 10 stress 0.1723312
## Run 11 stress 0.1754331
## Run 12 stress 0.1843778
## Run 13 stress 0.1705571
## ... New best solution
## ... Procrustes: rmse 0.0004753476 max resid 0.003322762
## ... Similar to previous best
## Run 14 stress 0.1842893
## Run 15 stress 0.172415
## Run 16 stress 0.1790648
## Run 17 stress 0.1864631
## Run 18 stress 0.1763068
## Run 19 stress 0.1778902
## Run 20 stress 0.1879778
## *** Best solution repeated 1 times
permanova(ps.sub, 'Household + Sex', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 11.7188 0.68925 2.2419 9.999e-05 ***
## Sex 1 0.0561 0.00330 0.5156 0.9684
## Residual 48 5.2272 0.30745
## Total 97 17.0021 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.sub, 'Household + Sex', 'wunifrac')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 4.8365 0.67991 2.1305 9.999e-05 ***
## Sex 1 0.0068 0.00096 0.1443 0.9985
## Residual 48 2.2701 0.31913
## Total 97 7.1135 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permdisp(ps.sub, 'Sex', 'bray')
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.01527 0.015267 0.9824 10000 0.3235
## Residuals 96 1.49198 0.015541
permdisp(ps.sub, 'Sex', 'wunifrac')
##
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 10000
##
## Response: Distances
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.00068 0.000684 0.0672 10000 0.8027
## Residuals 96 0.97753 0.010183
ggplot(sam) +
geom_line(aes(x = as.numeric(Household), y = as.numeric(Age..months.), group = Household)) +
geom_point(aes(x = as.numeric(Household), y = as.numeric(Age..months.), colour = Epileptic.or.Control)) +
xlab('Household') + ylab('Age')
permanova(ps.sub, 'Household + Age..months.', 'bray')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 11.7188 0.68925 2.3061 9.999e-05 ***
## Age..months. 19 2.1074 0.12395 1.0477 0.3502
## Residual 30 3.1760 0.18680
## Total 97 17.0021 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
permanova(ps.sub, 'Household + Age..months.', 'wunifrac')
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 10000
##
## adonis2(formula = model, data = df, permutations = permutations, parallel = core, strata = strata)
## Df SumOfSqs R2 F Pr(>F)
## Household 48 4.8365 0.67991 2.2162 9.999e-05 ***
## Age..months. 19 0.9130 0.12835 1.0569 0.3676
## Residual 30 1.3639 0.19174
## Total 97 7.1135 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1